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PREFACE 


Mathematical methods for computing satellite orbits consist of 
an algorithm of formulas - that can be programed on an automatic com- 
puter. This report is concerned, specifically, with methods for 
predicting satellite orbits; that is, given the satellite's state 
vector (usually position and velocity) at an epoch, compute the 
state vector at another epoch (or sequence of epochs). A large 
variS'ty of methods and techniques are available today for solving 
this problem and it is necessary to choose one that^ best suits a 
particular application. To do this, the methods must be converted 
intp. an executable program code and then analyzed to determine their 
efficiency and accuracy for producing usable results 1 

Each method consists of three basic parts: 

a. Formulation of the differential equations 'of motion, 

b. Analytical o'r numerical method for the solution 
(integration) of these equations, 

G. Mathematical model of the perturbing forces acting on the 
satellite . 

Each part can be coded as a separate module or subroutine in a com- 
puter program. The formulations (part. a) are discussed i-n section 
1.0 of this report-- New formulations are available that have a 
profound impact on increasing the efficiency of orbit predictions. 
This is due to the new stabilized differential equations that are 
easier to solve. The methods of integration (part b) -are also dis- 
cussed in section 1.0. The aim is to match formulation and inte- 
gration method for a particular application. Mathematical force 
models (part c) are discussed in section 2.0. 

In section 1.0 of this report, the concern is with “numerical 
accuracy" in the orbit prediction- This is defined to mean the 
accuracy to which the differential equations are solved, assuming 
that the programed equations perfectly describe the physical 
situation. The most recent advances in celestial mechanics and 
numerical analysis are examined. A set of rigorous numerical 
experiments were carried out to test the various combinations of 
formulation and integrator. 

The -subject of force models is concerned with "physical ac- 
curacy" of an orbit predictor. That is, how well does the predict- 
ed orbit agree to physical reality? This depends, usually, on 
how well the forces can be modeled. Several types of perturbing 
forces that affect the orbit are discussed in section 2.0. 

Questions to be considered are: 
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• What are the typical magnitudes of each perturbation? If 
they are neglected or poorly modeled, how is the accuracy 
of an orbit prediction affected? 

• ' How can the force, routine's be made less costly in 'terras of 

computer runtime and storage? What new,_ more efficient 
force models should be used? 

e What are the necessary force .models for a variety of 
standard Earth satellite orbits? 

The suggestions and recommendations . contained in this report 
are based on the particular cases, studied. An attempt was made to 
consider a variety of orbits, prediction intervals, output require- 
ments, etc. In view of the wide, variety of satellite missions that 
are being' planned, the scope of this study is necessarily limited. 
Nevertheless, it is hoped that these results and conclusions will 
be a useful guide in choosing the components of ah orbit predic-tion 
program, as well as a help in analyzing the numerical results of 
such a program. 


iv 



COMTENTS 


S.ection Page 

1.0 FORMULATIONS AND INTEGRATION METHODS 1 

1.1 Introduction ....... 1 

1.2 Formulations . 2 

1.2.1 Coordinate formulations ^ 

1.2.2 Variation of parameters ■ 5 

1.2.3 Total energy elements 6 

1.2.4 Choosing a formulation 7 

1.3 Numerical Integrators ^ ^ 

1.3.1 Adams formula 8 

1.3.2 Runge-Kutta formulas 9 

1.3. 2.1 Classical fourth~order 

formula .... - 10 

1.3. 2. 2 Fehiherg's formulas', ..... 11 

1.. 3.2.3 Bettis* improved formula ... 11 

1.3. 2. . 4 Shank’s formula 11 

1.3.3 General remarks 11 

1.4 Formulation-Integrator Combinations 12 

1.4.1 ■ Formula'tions investigated' 12 

1.4.2 Numerical integration methods 

investigated 13 

1.4.3 Comparisons ■ based on' storage and 

' cycle time 13 

1.4.4 Comparisons based on accuracy and 

execution time . ■ 14 

1.4. 4,1 “’Near-Earth o'rbit 14 

'1.4.4.2 Geosynchronous orbit 15 

1.4. 4. 3 Elliptical transfer orbit ... 15 

1.4.4.4 Highly eccentric orbit .... 16 

1.4.5 Conclusions from compari4ons 16 

1 .5 Analytical Solution Methods 21 

1.6 Kultirevolution Methods 22 

1.7 Concluding Remarks 25 


V 



Section Page 

2.0 MATHEMATICAL MODELS OF THE PERTURBING FORCES . . , 26 

2.1 Introduction 

2.2 Atmospheric Density Models ' 28 

2.2.1 Description of the models , 29 

2.2.2 Method of evaluation.. 31 

2.2.3 Effects of the Sun on the upper 

atmosphere 32 

2.2.’4 A new analytical atmospheric density 

model 33 

2.2.5 Orbit prediction , experiments ...... 37 

2.2.6’ Additional remarks ' 

2.3 Geopotential Model . . . ’ 1 ‘ . . ' . . . 'll 

2 . 3.1 Numerical exp.eriments . . I 42 

2 . 3.2 Theoretical background-. 44 

2. 3 . . 3 Comments and recommendations 45 

2.4 Luni-Solar- Gravity Mo.del’s 45 

2.4.. 1 Magnitude .o.f l, uni-solar ' 

perturbations 48 

2.4.2 Analytical model v.e.rsus- JPL .stored 

data 50 

2.. 4.-.3. 'Additional, sugges.tions • , • 57 

2.5 Time and Coordinate System Models 57 

2 . 5.1 Disoussion-. 58 

,2.. 5, ..2- Numerical exp.eriments 60 

2 . 5.3 Conclusions 63 

.2.6 Recommendations Based oh, Orbit Types 63 

2.6.1 Near-.Earth. .orbits. .64 

2.6.2 N.e.ar.-Earth orbit .lifetime studies ... 6'4 

.-2. 6,.-3 .Elliptical transfer orbits 65 

2.6.4 Geosynchronous orbit’s . 65 

2.6.5 General results 66 

REFERENCES 68 


v,i 



TABLES 


Table Page 

I KSMULT TEST RUNS 

(a) Prediction interval = 125 days . 23 

(b) Prediction interval = 165 days 24 

II DENSITY OUTPUT, JACCHIA VERSUS USSR 31 

III POSITION DIFFERENCE, JACCHIA VERSUS USSR 32 

IV POSITION DIFFERENCE, JACCHIA VERSUS AMDB, 

NODRAG 33 

V ORBITS USED IN PREDICTION EXPERIMENTS 37 

VI POSITION DEPENDENCE ON DENSITY MODEL 

(a) Orbit A 38 

(b) Orbit B 38 

(c) Orbit C 39 

(d) Orbit D 39 

(e) Orbit E 39 

(f) Orbit F 40 

VII EXECUTION TIME AND STORAGE REQUIREMENTS 4l 

VIII POSITION DEPENDENCE ON GEOPOTENTIAL MODEL 

(a) 2x0 versus 2x2 43 

(b) 8x0 versus 8x8 .44 

IX ORBITS USED IN PREDICTION EXPERIMENTS ■ 47 

X LUNI-SOLAfl GRAVITY EFFECTS ON POSITION 

(a) Combined Sun and Moon ............ 48 

(b) Sun only 49 

(c) Moon only 49 

XI POSITION DIFFERENCE, JPL VERSUS ANAL 50 

XII COMPUTER RUNTIME AND STORAGE, JPL VERSUS ANAL ... 57 


'Vii 



Table 


Page 


XIII 


XIV 


TIME AND COORDINATE SYSTEM COMPARISONS 

(a) GNP versus GNP* 

(b) GNP* versus GNP** , . . . . . . 62 

EXECUTION TIME COMPARISONS (GNP) 63 


viii 



FIGURES 


Figure 

1 

2 


3 

4 

5 

6 


Independent variables . 

Numerical efficiency curves 

(a) Near-Earth orbit 

(b) Geosynchronous orbit . . . 

(c) Elliptic transfer orbit 

(d) Highly eccentric orbit . . . 

Atmospheric density variations 

(a) Variations with altitude and epoch 

(b) Variations with position of the Sun at 

500 km altitude 

Position of the Sun, JPL versus ANAL 

(a) Right ascension ' . . . 

(b) Declination 

(c) Radial distance ... 

Position of the Moon, JPL versus ANAL 

(a) Right ascension ...... 

(b) Declination . . . . 

(c) Radial distance ..... 

Precession nutation, and Greenwich rotation . . . 


Page 

3 


17 
1 8 

19 

20 


34 

35 


51 

52 

53 


54 

55 

56 

59 



APPLICATION AND ANALYSIS OF SATELLITE 
ORBIT PREDICTION TECHNIQUES 
By Analytical and Computational Mathematics, Inc. 

1.0 FORMULATIONS AND INTEGRATION METHODS 


Introduction 

The computation' and prediction of sa*t.ellite- orbits requires' 
the solution of a set of ordinary differential equations. Studies 
have been carried out to determine the best methods for obtaining 
these solutions. It will be assumed in this section -'that the per- 
turbing force's are known exactly. The goal, therefore, is to in- 
vestigate the numerical accuracy of a satellite orbit' computation 
program. That is, to determine the effects of roundoff and 
truncation errors on the solution. 

Many different formulations of the satellite differential equa- 
tions are available. Of speciai interest here are the new formula- 
tions that have the mean motion based on the total energy. They 
will be compared to the more classical formulations and evaluated 
via numerical experiments.' 

Me-thods' for the numerica'l solution of ordinary differential 
equations are also discussed. Numerical experiments have been 
carried out to determine the most efficient methods, in terms of 
computer cpu times. 

An intention of these studies has been. to determine the formu- 
lation and numerical integration combinations which exhibit most or 
all of the following attributes. 

a. Numerically Accurate : The combination should produce so- 

lutions which are consistently more accurate than the force model 
accurac.y (see section 2.0 on force model errors). In other words, 
roundoff and trunca'ti.on errors should be much smaller'* than modeling 
errors. 

b. Predictable' Accuracy : In .addition to point (a), the accu- 

racy of, the solution should be -predictable . In practical applica- 
tions i, it is nece'ssary, that an uppe'r bound on the error be known. 

c. Efficient ; For a given accuracy,' a combination should 
require as little computation time as possible to' obtain the 



solution. Thus- i‘t is desirable to have a method which is 
reasonably accurate but requires the least computation time. 

d. Reliable ; A combination should give reliable results for 
all possible applications. In addition, the combination should 
have a "built-in" reliability so that it is insensitive to possible 
misuse of the algorithm. 

e. Stable I Small errors in the initial state should not 
cause the combination to run into numerical difficulties. Some 
combinations may be reliable over short propagation times; but, 
because of instabilities become completely erroneous for longer 
times. 

f. General ; The formulation should be valid over the whole 
range of intitlal conditions. For instance, certain sets have 
singularities for circular or equatorial orbits. However, since we 
are most int-erested, in elliptical orbits, we wi.ll exclud.e the gen- 
erality to parabolic or hyperbolic orbits. 

S' Concise : The computer algorithm for a combination sho.uld 

require a small percentage of the storage requirements of the total 
force model , 

A brief discussion is given on analytical and seminumerical 
methods. These are specialized techniques that allow a large sav- 
ings in computer time. Their accuracy and r^ange. of applications 
willbediscussed.- 


1.2 Formulations 

The satellite differential equations can be .formulated in a 
number of different ways for different purposes. A‘ll formulations 
fall basically into three classes: 

a. Coordinate Formulation 

b. Variation of Parameters 

c. Total Energy Elements 

These three .classes can each be divided into' two' sets, formula- 
tions with time as the independent variable and formulations 
with an independent variable other than time.' Before discussing 
the different classes one should understand the effect of the 
independent variable on the numerical integration. The ellipses' 
of figure 1 show where the force model is_ being evaluated for 
equal steps in three different independent variables. "The perigee 
point, of this ellipse (e = .7) is located on the right-hand side. 
Note that equal steps in the true anomaly tends to evaluate the 
.force model more near the perigee point. Onthe other hand, 
equal steps in time tends to bunch the majority of the /force 
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True anomaly 




.Figure 1.— Independent variables. 
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evaluations near the apogee. The eccentric anomaly smooths 
the evaluations equally over perigee and apogee. Since the 
stronger perturbations are usually associated with the perigee 
point (atmospheric drag, for example), it becomes immediately' 
clear why the use of the true or eccentric anomaly results in 
an automatic "analytical” step size control. 

One could conclude from figure 1 that if the perturbations 
are much stronger at perigee than at. apogee, then the true anom- 
aly would be the most appropriate choice for independent variable. 
If the perturbations at apogee' and, perigee were nearly equal then 
one might expect the eccentric anomaly to be the better choice. 
Clearly, time would be a poor choice for either case. Of course, 
when the orbit is circular then the different independent variables 
become quite the same. 

1.2.1 Coordinate f nrmii lati on ■ - Let the coordinates x 
of the satellite be referred to an inertial rectangular coordinate 
system. The most basic formulation of the equations of motion 
is the so called Cowel method (ref. 1) 


V 3V . ^ 

+ p 

r3 2 k 


where 


V is the perturbing potential function and P represents acceler- 
ations that are not derivable from a potential. The gravitational 
parameter is p. This method uses time as the independent vari- 
able. The main advantage is, of course, it's simplicity. The 
approach, however, has several disadvantages. The right-hand side 
of the differental equations' are quite large- and therefore one is 
forced to take small numerical steps or use a high order integrator. 
Also, the differential equations are unstable. An attempt to alle- 
viate the first problem was made by Encke (ref.,1). If the equation 
of unperturbed motion is 


-»• 

t + y — = 0 

and 

Ax >= X - ? 
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Then the Encke equations read 



The right-hand side of the equation is small if Ax is small. But 

•¥ 

Ax may grow to be quite large, therefore, this method must be '*recti- 

fie.d." This rectification involves setting C = ’x and &x = 0. This 
oan cause additional instabilities in the differential equaitions. 

One .may modify the Cowell 'equations by using a fictitious time 
s,. defined by the differential equation 



The equations of motion then become 



This method (ref. 2) retains the coordinate formulation and adds an 
important analytic stepsize control, but it requires an extra dif- 
ferential equation to compute the time. An Encke-type formulation 
can also be applied to these equations to reduce the magnitude of 
the right hand side of the differential equations. 

1.2.2 Variation of parameters .- By solving the unperturbed 

problem one obtains a vector of six constants of integration a- 
Classical variation of parameters (VOP) methods develop equations 
of motion which define how the six constants of integration vary 
in the perturbed case. They have the form 

a = 

The r.ightrband sides of the differential equations are small, as 
in the Encke method, but do not require the expense of rectification. 
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The Lagrange and Deluanay element formulations are typical VOP 
examples. These sets, however, have singularities at zero incli- 
nation and zero eccentricity. Another example, the Poincare' 
elements, avoids these problems. These VOP elements, which prove 
very fruitful in analytical studies, are somewhat cumbersome for 
numerical Integration. Since the forces are usually given in 
Cartesian coordinates, the differential equations require a con- 

version from a and t to x. Another method, developed by Pines 
(ref. 3), simplifies this conversion. This, method, developed with 
time as the independent varianie, has mixed secular terms. These 
terms grow in magnitude with time and will eventually degrade the 
accuracy of the solution. The secular terms were later eliminated 
by a conversion to a perturbed time as independent variable (ref. 4). 
Another disadvantage of these classical variation of parameter meth- 
ods is that a costly iteration to solve Kepler's equation is required 
at each integration step (ref. 5). 

Burdet has developed a VOP method based on Sperling's solution 
of a set of regular and linear differential equations of the unper- 
turbed (two-body) motion (ref. 6). Although this approach requires 
14 differential equations, it provides increased numerical stability. 
The method presents no mixed secular terms and requires no Kepler 
iteration. The regularization requires a transformation to an inde- 
pendent variable other than 'time. For elliptical cases the method 
can be specialized so that the independent variable becomes the 
eccentric anomaly and thus this method also has an analytical step 
size control (ref. 7). 

1.2.3 Total energy elements .- Kustaanheirao and Stiefel have 
also developed a linear and regular set of differential equations. 

The resultant differential equation is that of a perturbed harmonic 
oscillator in 4-space 

P' + m ^ 


where 



However, the major difference between the Burdet and KS method is 
that the total energy, instead of the two-body energy, is used as 
a parameter. This has been shown' by' Stiefel and Scheifele (ref. 8) 
to significantly improve the accuracy and stability of the solutions. 
The KS method, which has ten differential equations, can be special- 
ized for elliptical orbits in which ease the eccentric anomaly be- 
comes the independent- variable. Variation of parameters applied to 
the integration. constants of the four dimensional harmonic oscillator 
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provides eight element differential ^equations. _ Two addi- 
tional element equations come from the perturbed fr.equen.cy 
tbased on tne total energy) and an element t for the perturbed 
time . 

The Burdet elements have been modified by Bond (ref. 9) so 
tnat they, too contain the total energy element. Both element sets 
^KS and Bur dets-Bpnd ) have -proven to be extremely accurate, stable, 
and yet have concise formulations. The two methods, however, are 
'not well sui.ted, for analytical solutions. 

Recently, a number of canonical' element sets have been devel- 
oped (refs. 10, 11, and 12) which incorporate a total energy element 

in an extended phase space. iney nave eight eanonical dependent 
variables with an independent variable other than time. Two oi the 
sets, DS<|i and DSu, are similar to the Delaunay elements and have 
the,- true and eccentric .anomaly, respectively as their Independent 
variable,. These sets have singularities for small, eccentricities 
and inc li'ha-tions . They have been transformed, to Poincare' type 
.elements,- PS'}', -and PSu, to remove the s.ingularlti es . Unlike the 
classical, e.lements,, ,no iteration of Kepler's equation is required. 

All of -these "extended phase space" .sets have b.ee.n s^hown to be as 
accurate .and. efficient as the KS set,, .but are not as concise in 
their formulation.. Becaus’e they are canonical sets, they are more 
readily adopted to analytical perturbation .theory. 

1.,2„4 -Choosing -a formulation .- In examining the different ele- 
ment fqrm.ulat‘ioris , one finds a number of characteristics which are 
considered desirable f.or use in a numerical orbit prediction pro- 
gram. 

a. The elements should vary slowly and .smoothly as a function 
of the independent .variable . 

b. The elements ' should always be defined for elliptical orbits. 

c. The independent variable should be the true or eccentric 
anomaly, allowing for analytical step size control. 

d. The energy element should include the total energy. 

e. rne formulation should require, no iteration of Kepler's 
equation during the integration procedure. 

Only the KS, Burdet-Bond, PSif) and PSu elements of the total 
energy class exhibit all of these characteristics. Indeed, a wide 
range of experiments have consistently shown that the total energy 
formulations, when compared to the more classical formulations, are 
faster, far more stable, and yield more accurate results, while re- 
quiring very little (if any) more computer code (ref. 13 and 14). 

Two of these four sets, KS and PSifi , were chosen for further 
testing and comparisons. The results and conclusions from the nu- 
merical experiments are discussed in section 1,4. An additional 
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objective was to match a numerical integration method to these 
two formulations. Several numerical integrators are discussed in • 
the next section. 


1.3 Numerical Integrators 

Given the differential equations and initial conditions', one 
may employ finite calculus to numerically integrate the differen- 
tial equations, giving the solution as a function of the indepen- 
dent variable and initial conditions. Two classes of integration 
methods will be discussed in this section: 

a. Runge-Kutta (single-step) formulas 

b. Adams (raultlstep) formulas 

The major difference between these two classes is that the single 
step methods develop the solution over one step, requiring no pre- 
vious information on the solution. The multistep formulas require 
a sequence of values of the derivatives in order to advance the so- 
lution one step. If the derivatives vary slowly and smoothly from 
step to step, then both classes perform well. I-f the derivatives 
become erratic or impulsive, the single step methods, with step size 
control, are preferred. Step size eontro'l with the multistep method 
is generally not possible. 

Note that a step size control procedure that is built into an 
integrator is called ’’numerical" step size control. Compare this to 
the "analytical" step size control that is introduced in section 1 . 2 , 

1.3.1 Adams formula .- The classical Adams method is an example 
o.f a multistep formula. It uses the Adams-Bashf orth formula (pre- 
dictor) and the Adams-Moplton formula (corrector) to obtain the solu- 
tion. A complete discussion of this method is given by Henrici (ref. 
15 ) . 

The solution is extrapolated one step forward by applying the 
predictor formula. The derivative is then evaluated with the pre- 
dicted solution and a corrected solution is obtained from the cor- 
rector formula. This classical pr edict-correct type of algorithm 
is not well suited for variable-step applications. Therefore, only 
fixed-step Adams integrators are included in this study. A trun- 
cation error estimate is available and can be used to determine the 
best step size. 

Consider a vector that is defined by the first-order vector 
differential equation 
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where x is the independent variable. Also define 
X. = X + jh 

JO 

yj - y(xj) 


£3 = f(.j.yj) 

The predictor and corrector formulas are 

N . 

n+1 '■'n k n4-k-N 

k=0 
N 


y. 


n+1 


" ^n+k- 


■N+1 


k=0 


where N+1 is the order of both formulas. 


Observe that th,e corrector formula is ah implicit equation 

since ^^+1 depends on Y^+l * Usually, one iteration is sufficient 

for solving the equation. ' The coefficients are functions 

of the order of the formula. These coefficients may be compute.d re- 
cursively for any given order before the integration is initiated. 

For this reason the Adams method is sometimes called a "variable order" 
method . 

The truncation error of the corrected solution is given by 



where is a function of only the order. 

The first N + 1 values of f / are required to start the 

integration. In the numerical examples in Section 1.4, an eighth 
order RK formula was used to start the Adams integration. 

1 . 3.2 Runge-Kiitta formulas .- Runge-Kutta (RK) formulas fall 
in the class of single-step integration formulas. That is, they 
are self-starting. In general, an RK formula can be written as 

V 

y =£ y +y^ W. K. 

^n +1 ^ 1 1 

i=l 
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where are weighting coefficients, v is the number of func- 


tion evaluations required at each step and is define,d as 

“ (*n + 'i'*’ y,. 

\ J-1 

Equations of condition for W^, and A^j are found by 




expanding the RK formula in a Taylor series. Solution of these 
equations provides the numerical values of the coefficients. To, 
date, however, there is no simple manner to obtain values of these 
coefficients for any order. A good discussion on RK methods is 
contained in reference 16. 


A short description of five different RK formulas is given 
in the following subsections. A list is given" below showing the 
number of "function calls" (NFC) that are required at each step 
by the various RK formulas. A function call is one evaluation of 
the differential equations. 


RK ■.Co.r.qmXa 


NFC per step 


RK4 

RKF45 

RK45 

RKSS 

RK78 


4 

6 

6 

10 

13 


1.3. 2.1 Classical fourth-order formula; The classical fourth- 
order formula (RK4) is 


'n+l 




+ 2K^ + 2K^ + 


.) 




Kj - h£ 


Kj-hf 




X -J- i'h, y + i K . , 
n 2 n Z Z I 


K. 


== hf(x^ + h, + K3) 
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It is' 'very 'po'pular and quite efficient. However, when compared to 
other available RK formulas, it has the following disadvantages. 

'a. No step size control. 

b. Less efficient for accurate integration. 

Its chief advantage is the small computer_ storage requirement. 

1.3-2.2 Fehlber.g's formulas: Fehlberg's RK formulas overcome 

the two disadvantages mentioned above. His formulas contain an auto- 
matic truncation error estimate and are available in orders one 
through eight (refs. 17 and 18). Reference 19 contains a discussion 
of the error propagation, as well as the coefficients for each RK- 
Fehlberg formula. A fifth order version (RKF45) and an eighth order 
version (RK78) have been tested. 

r 

1 ..3 . 2 . 3 Bettis' improved formula: New versions of Fehlberg's 

formulas are being developed by Dale Bettis at the University of 
Texas, Austin. He uses a Davidon optimization scheme to compute new 
RK co'ef ficients that minimize the truncation error. The expression 
for the truncation error is, known as an ana'lytical function of the 
coefficients. Therefore, the optimized coefficients apply to any 
set of differential equations and, in effect, raise the order of the 
formula (ref. 20). 

Dr. Bettis has made available to the. Mission Planning and Anal- 
ysis Division .(MPAD) his optimal version of RKF45, which will be 
referred to here as RK45. There is no increase in computer run time 
or storage over that for RKf 45 since the algorithm is the same. 

1.3. 2. 4 Shanks formula: An additional RK formula developed by 

Shanks (ref. 21) has also been tested. It was chosen since it was 
of high order and has been shown- by other* investigators (ref. 22) to 
have accuracy comparable to that of Fehlberg's formulas. 

Shanks attempted to minimize the number of function evalua- 
tiions per step for a given order. The formula tested here (RKS8) 
is approximately an eighth-order formula and requires only 10 func- 
tion evaluations per step (RKF78 requires 13) but it contains no 
step size control. Therefore, it is expected that RKS8 will be 
more efficient than RKF78 when a constant step is desirable, such 
as for circular orbits. The reverse situation is expected for very 
elliptical orbits. 

1.3.3 General remarks '.- For the special case of orbital motion, 
the right-hand side of the differential equations are usually smooth 
functions of the dependent variable. In this regard, the Adams for- 
mulas are expected to be more efficient • than the Runge-Kutta formulas. 
Elliptical orbits may,' however, require the variable step size option 
of the RK formulas. 
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The required order of an integrator for a given appuracy de- 
pends on the formulation of the.' dif f"erential. equations. A method' 
such as the Cowell formulation requires a high order integrator to 
obtain satisfactory results. But, a method suph as the KS formu- 
lation is somewhat insensitive to the order of the integrator (see 
section I.M), i.e., all integrators perform about the same. 

Altho.ugh high order integrators are usually more accurate 
they are also less stable (ref. 22). Thus, in choosing the inte- 
grator. one must find a happy medium where the order is large enough 
for reasonable accuracies but small enough not to introduce instabil- 
ities. 

1,^ Formulation-Integrator Combinations 

When building a numerical orbit prediction program it is de- 
sirable to choose the formulation and integration method that best 
satisy the requirements discussed in section T. 1 (Force model require- 
ments are discussed in section 2.0). The discussion in this section 
concerns the choice of these two important components and relies on 
both theoretical considerations and numerical experiment's. Use was', 
made of the results and experience of prior investigators to eliminate 
some methods and limit the number o'f possible combinations.' 

Several different problems of orbit, prediction were chosen for 
the numerical experiments. The orbits chosen are important for 
shuttle orbiter and payload missions. Also, they demonstrate the 
different qualities of the 'dif f erential”equation formulation -- numer- 
ical integration (DE - NI) combinations.' A thorough discussion of 
the comparison is contained in reference 23. The following gives- a 
summary and the conclusions. 

1.-4.1 Formulations investigated .- 'The discussion in section 1.2 
concluded that there are four available element sets (KS, Burdet-Bond, 
PSff, PSu) that have significant advantages over' all the rest. These 
four are in the class of "total energy elements" and each. is 
based on an independent variable different from time. It was de- 
cided to choose two of these for further testing and evaluation. 

In choosing between the four total energy element sets, se.v- 
eral considerations were. made. First, PS^ has true anomaly as 
the independent variable, the remaining three have eccentric' anom- 
aly. Thus, PS(j> was chosen for further testing in order to deter- 
mine the advantage of the step spacing (fig. 1) of the true anomaly. 
For orbits that are strongly perturbed by the geopotential and drag, 
PSiJ) should -hav.e an advantage. 

KS was chosen from among. PSu and .Burdet-Bond because (com- 
pared to PSu) it has a concise formulation and it has fewer ele- 
ments than Burdet-Bond (10 versus 13). Numerical experiments and 
theoretical analyses of Bond (refs. 10 and 24) indicate that the 
PSu has a slight advantage in long term stability. However, this 
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does not seem to be significant enough (for most applications of 
orbit predictors) to offset the value of a concise formulation. 

1.4.2 ' Numerical inteRration methods investigated .- Graf 
. 13) investigated a number of different numerical integrators 
the Cowell and KS formulations-. From the conclusions of that 
rt it was determined that four integrators had advantages 
warranted a thorough investigation with the KS and PS(t) formu- 
ons . A short descrip'tron of these integrators and the reasons 
they were c,hosen are listed below: 

a. AD9 - This is an' Adams (multi'step) method of ninth ord 
which evaluates the deriva'tive once per step. This- fixed step 
od proved to be very efficient with the KS method even for eoce 
orbits. 

b. RK4 - This is the classical four.th order 'Runge-Kutta 
(single-step) method which requires four derivative evaluations per 
step. Even though this is a simple, low order method it proved effi- 
cient with the KS method. 

c. RK45 - This is the fifth order Runge-Kutta method with the 
option of step size' co.ntrol . It requires six derivative evaluations 
per step and uses Bettis's optimized coefficients because <of their 
proven accuracy (ref. 20). This method was chosen mainly because 

of i-ts powerful step size control. 

'd.. RK7'8 - This is the eighth order Runge-Kutta-Fehlberg meth- 

od that requires 13 derivative evaluations per step. It does have 
the option of step size control. It was chosen for its efficiency 
for stringent accuracy requirements. 


er 

meth- 

nt-ric 


(ref 

with 

repo 

that 

lati 

why 


1.4.3 ' Comparisons based Qn_^storage and cvcl'e time .- In com- 
paring '*KS and PS<S> '^one might argue that t-he KS formulation has a 
clear advantage because l.t Is so concise. Certainly, if one com- 
pares storage one finds that PS4i equations require 100 percent 
more computer storage than KS equa'tions. However, this is not 
qui.te a fair comparison. One should always compare- formulations 
by the percentage of the formulation storage to the storage re- 
quired by the total orbit prediction algorithm which includes the 
force model and numerical integrator. Thus with a force model 
which adequately represents the forces on a satellite about the 
Earth and with an integrator such AD9, KS composes 13 percent of 
the total storage and PS<I> only 23 percent. Therefore, the im- 
pact of the PS^ additional overhead is not as large as it might 
first appear, but still should be considered for very str.lngent 
storage requirements. 


Vie define cycle time as the time required for one evaluation 
of the- derivatives of an element set. Certainly KS will have a 
smaller cycle time than PSit). -But for accurate predictions, a 
large force model is required arid thus "most of the time evaluating 
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the derivatives is spent In. the force model. Jhus one finds that, the 
cycle time for ■ PS(j> is only 6. percent more than the KS cycle time.. 

One- can only conclude that, "for accurate orbit . predictions 
where, a large force, model is required, 'the overhead, of ‘ PSij) is 
negligible except for very stringent storage .and cycle time re- 
quirements - 


,1.4.. Comparisons based on accuracy and execution time .- The 
four integration methods and the KS and PS<> formulations' have been 
programed in double precl'sion on the Univac 1110 (ref. 25 -)-. The 
results. of the various, comparisons .are displayed in graphical form 
(,figs. 2,(a) through 2(d)). The bbt.tom scale in .each graph is the 
execution time (cpu time) i-n s.e.c.onds and 'the .left-hand scale is the 
accuracy defined as 


6 = -Iqgio^R ^ 


where 

aR is the position vector difference magnitude ,betw,een the 
test and reference solutions at the final time. . 

R is the positi,on magnitude 'of. the reference solution which 
is -aceurate to 10 digits. 

No.te .that, the cpu time was varied by. taking suoce.ssively smaller- 
step sizes .or. tolerance . criteria in each of the NI - DE combina'tions . 
Methods whose curves lie above and to the left are therefore' con- 
sidered the most efficient. 

Again, one should refer to reference 23 for a more detailed 
description of these numerical comparisofas, .graphical displays, and 
conclusions. 

1.4.4..1 Hear Earth orbit:- -This test .case has/ the foll.o.wlng 
init-lal -conditions-' 

- = 278 km- I = 28.4°. Epoch of Jan. 1, 1975 

hp = 167 0) = n = M = o'* ’ ' 

The force model includes air drag and an eighth order, eighth 
degree, geopotentiai model.. The final time of integration is 
= ,2.0 days which is about 32 revolutions. 


The KS and PS performed about the same, which could be ex- 
pected for this near circular orbit. All the Runge-Kutta methods 
performed about the same and ■ both . formulations appear, bo be in- 
sensitive to the order of the integrator. In both formulations. 
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the AD9 method was slightly more efficient than the FK methods. 
This too could be expected since the forces vary slowly and are 
very smooth in this case. Figure 2(a). displays the efficiency 
curve of the KS and PS combinations with AD9 and RKM5. 

Geosynchronous orbit: This test case has the fol- 

lowing initial conditions 

h, = h„ = 35862 km Epoch of Jan. 1, 1975 

a P 

f>=n=M = 0.0° 1 = 0,0° 

t^ = 100 days = 100 revolutions 

The force model includes the Sun and Moon perturbations and a 
fourth order, fourth degree geopotential model. 


Results ; 


Again, the KS and PS methods performed almost the same. For 
accuracies of or about 100 to 1000 meters, the RKll , RK45 

and AD9 are almost equivalent. For slightly more stringent accura- 
cies-, the AD9 integrator is most efficient in either formulation. 
The RK78 fails to be competitive except for very stringent accura- 
cies of 6>7 (.less than 1 meter). Figure 2(b) displays the curves 
of the KS and PS combinations with AD9 and RK4 . 


1 . 4 . h . 3 Elliptical 
this case are 


transfer orbit; 


The initial .conditions of 


hp = 200 km 
h., = 35862 km 

a 


I = 30.0° 

oj = = M = 0.0° Epoch of Jan. 1, 1975 


= 2,.0 days ~4 revs 


The force 
potential 
forces at 


model includes drag, an eighth order, eighth degree geo- 
and Sun-Moon perturbations. Note that in this case the 
perigee are much stronger than those at apogee. 


Results ; 


In this case the PSti> formulation 
efficient than KS, regardless, of Integra 
pected since PS4> uses the true anomaly 
able. Fixed step AD.9 combination with P 
ciency. But the SK45 with its excellent 
used with either KS or PS4* for competi 
RK78 showed. to be inefficient again exce 
Also the fixed step methods. .AD9 and RK4 
the KS formula-tion . This is because the 
(eccentric anomaly) is not well suited f 
displays the KS and PS combinations with 


proved to be much more 
tion method. ‘ This was ex- 
as the independent vari- 
S«i> showed the highest effl- 
step size control could be 
tive efficiencies. 'The 
pt for stringent aocuracies. 
fared poorly when used with 
KS independent variable 
or this case. . Figure 2(c) 
AD9 and RK45.'‘ Note the 



large difference between the AD9 combinations but small difference 
between, the RK45 combinations. 

1.4. 4. 4 Highly eccentric orbit: The initial conditions of 

this case are 

hp = 425 km I = 30® . 

h, = 258,903 km w = s M = 0.0° Epoch of Jan. 1, 1975 

Si 

e = 0.95 

The force model includes drag, an I8th order geopotential and the 
Sun-Moon perturbations. The final tLme of integration is t^ = 

50 days ~ 8.6 revolutions. Note in this case the perturbation 
forces are equally strong at perigee and at apogee. 


The KS formulation showed bo be the stronger formulation in 
this case which is again to be expected since the eccentric anomaly 
as independent variable is better suited for this highly eccentric 
orbit. The RK45-KS combination was the most efficient but one could 
use the RK45-PS and even the fixed step AD9-KS combination without 
much loss in efficiency. The RK78 integrator appeared to have sta- 
bility problems and its step size control did not perform well. 

Since the tru’e anomaly is not well suited with this case, the fixed 
step methods AD9 and RK4 did not do well with PS<j>. Figure 2(d) 
shows the accuracy curves of the PS<|> and KS formulations with 
RK45 and AD9. 


.1.4,5 Conclusions from comparisons .- It was concluded from 
previous experiments that the PS4> and KS formulations possess most 
or all of the attributes listed in the introduction. The comparisons 
have shown the KS method does have a slight advantage iji relation to 
storage requirements and cycle time. Results from the numerical 
studies show that both methods' are equally powerful for circular 
orbit oases but differences between the two formulations become 
apparent for more eccentric orbits. 


It is recommended for circular orbits with slowly varying 
forces (such as the two cases that were examined) that the AD9 
integrator be used with either the KS or PS4> formulation for the 
most efficient results. However, the RK4 and RK45 integrators may 
also be used with little loss in efficiency. If the force model 
includes discontinuous forces such as venting, one may find that 
a low order single-step method such as RK4 to be more adequate. 


In the case of eccentric orbits where the perturbing forces 
are much stronger at perigee than at apogee (such as the elliptic 
transfer orbit), the PS«j’ formulation should be used. For best 
results the AD9 integrator is recommmended . However, an integrator 
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(a) Near-Earth orbit. 

Figure 2.- Numerical efficiency curves. 
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(b) Geosynchronous orbit, 
Figure 2.- Continued. 
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(c) Elliptic transfee ofbit. 
Figure 2.- Continued,, 
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(d) Highly eccentric orbit, . 
Figure 2 .- Concluded. 
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with - good step' size control -such as RK^I’5 may be used with KS’ 
formulation for competitive- result's*. ' 

For highly eccentric orbits, where the perturbing forces are 
important at both perigee and apogee, the KS formulation is best. 
Although the KS formulation has an analytic step size control that 
is suited for this case, the method produces the most efficient 
results w'ith ’the' help of th'e excellent RKM5 numerical ■ step size 
control. Ad'equate results may also be’, obtained from the RK45-PS 
combination. 

From the s.tudy, it appears there is no clear cut "winner" 
between KS and PS4> formulations. Certainly, both methods are 
very’ powerful and should be considered 'as the’ leading contenders 
for numerical orbit prediction algorithms. The AD9 method cer- 
tainly is the best integrator when the analytical step size con- 
trol is adequate but the RK45- becomes powerful when additional 
help is needed for step size' control . 


1.5 Analytical ’Solution Methods 

Whereas the numerical methods that were discussed in the pre- 
vious sections develop the solution in discrete steps, analytical 
methods produce the solution in the form of finite mathematical 
expressions. These ' expressions contai'n explicitly the dependent 
variable. Numerical- evaluation of these’ expressions provides the 
position and velocity of the satellite at a given time- 

An analytical solution can be thought of as a ".one step" 
method. When used for orbit prediction,' the ‘computation cost to 
obtain the state is always the same, regardless of the prediction 
interval. ’ T*ypical computation time is less than .one second.’ 

Analytical solutions to the J 2 (oblateness) problem have 

been established In the DS-elenients (-ref. 11) and PS-el-ements " 
(ref. 12). An orbit prediction program (ANALYT-)- based on the DS- 
ana'lytical solution has been developed and is documented in 
reference 26. Numerical" experiments contained in reference 26 
show, that ANALYT solutions of the ^2 problem have errors on the 

order of a few meters, so long as the ecoentri'city is larger than 
0.01. In addition, this error remains constant over several hun- 
dred revolutions (ref. 11). 

.> 

Additional testing -of the ANALYT program was carried out in 
reference 27. The refer'ence tra jeotori'es ’were obtained from the 
KSFAST program (ref. 28), using an l8th' order,- 18th degree geo- 
potential model plus' atmospheric drag. It was found that the pro- 
gram gives good results for those orbits where drag is not so 
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important; i.e., height of perigee is above 700 km. Position errors 
for a prediction of 60 revolutions was on the or.der of a few kilo- 
meters - 

1.6 Multirevolution Methods 

t 

A seminumerical orbit prediction method has been developed, 
based on the - mul'tir evp-lu tion integration technique.. The method 
makes use of the fact that the orbital motion of a satellite is. 
nearly periodic from revolution to revolution, as measured from 
some orbital reference point such as perigee.- The algorithm can 
extrapolate the satellite’s orbital- elements many revolutions 
ahead, thereby saving much computation time. Typical applications 
are for lifetime studies, orbit stability analyses and reference 
tra j ectori es . 

Let yp be the osculating elements at a prespecified orbital 

reference point. Then the multirevolution method solves the fol*i 
owing first order difference equation, 

7p+1 “.yp = f(yp.p) 

where ..p is the revolution number. The above equation results 
from the finite change in the elements over one revolution, meas- 
ured from the reference point. This change may be computed by a 
numerical integration of the equations of motion, with the use of 

yp as initial conditions. Since, p will have integer values only, 

j ^ ..... 

the elements <=y are known at discrete points of a continuous .inde- 
pendent variable, forming a grid of equal intervals. The multi- 
revolution method provides the solution of the difference equation 
at points separated by M grid points, where M is the multirevo- 
lution "stepsize". The computed solution, is then established on a - 
large grid, each interval of which contains M subintervals. 

The multirevolution formula is written as 

N 

^(n+l)M "" ^nM S “j ’ 

j=o . ■ ■ 

P = 0, 1 

The subscripted expressions on y 
number. p = 0 gives a predictor 
formula.- n is the step number, 
defined as 


and f refer to, the re.volution 
formula and p = 1 a corrector 
The- backward differences V* are 
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is very ^similar to the -Adams formula (section 1.3.1). Further 
theoretical details of the multirevol.ution algorithm are given in 
references 29 and 30 , ' 

An orbit prediction program KSMULT' based on the multirevolu- 
tion algorithm has been developed and is documented in reference 31 . 
The KS elements are used for' extrapolation and the differences are 
computed by the KSFAST program. Therefore, all of -the input/output 
and. force, .mod ei options of KSFAST are available to KSMULT.. 

Numerical evaluations of KSMULT hav-e been carried out (ref. 32). 
The purpose was to determine the optimum values for the order N 
and stepsize M, when the multirevolution algorithm is applie.d to 
ne.ar Earth orbits. The initial parameters of the orbit investigated 
were: 

altitude 296 km 

eccentricity 0 

inclination 3,0'°‘ 

period 90.^1 min 

Perturbations included an, 1.8th order zonal geopotential model plus 
air drag (Jacchia density .model). Re suits' are shown in table. 1(a) 
and (b) for prediction, intervals of 125, days and 165 days, respec- 
tivel-y. 



(a) 

TABLE I.^ 
Prediction 

KSMULT test; RUNS 
Interval = 125 days 


li_ 


AR- 

CDU time 

Savings 

6 

re 

2 ,.. 0 km 

162 .sec 

6. '8 

6 

24 

l 'l ..1 

118 

. '9.- 3 . 

8 

re 

.8 

169 

6.5 

10 

‘l6 

8.9' 

175 

6.3 
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TABLE I 


Concluded 


(b) Prediction Interval = 165 days 


Sa vings 


M 

v\R 

cpu time 

factor 

16 

9.1 km 

204 sec 

7.1 

16 

1 .6 

211 sec 

6.9 


The reference solutions were obtained from KSFAST . 
iR is the vector magnitude of the position difference 
in the KSMULT and KSFAST solutions. This is primarily 
a downrange error. 

The errors in these predictions are’ oh the order of a 
few kilometers. This is less than the expected errors 
due to force model inaccuracies (see Section 2.0). 

Ratios of KSFAST opu time ’are shown in the column under 
"Savings Factor." This gives an indication of the 
savings to be realized with KSMULT. Note that the total 
KSMOLT runtime for precision computations’ of 2628 revo- 
lutions (table 1(b)) is only about 3 1/2 minutes. 

According to these results (and addtional results in 
reference 32), M=16 and N=8 are the optimum values 

for near Earth orbits. 

This particular orbit decayed at 170 days. During the 
final stages of the decay, the orbital elements change 
rapidly and multirevolution is no longer the appropriate 
numerical integration method. The program can then 
switch to KSFAST for computation of the final few revo- 
lutions . 

An additional multirevolution program STEPR has been built, 
based on the routines in the general orbit prediction program (GOPP) 
(ref. 25). STEPR is described and documented in reference 33. All 
of the element sets (KS, PS<|), PSU, DS4>-, DSU , Cowell) and inte- 
gration methods (Adams, RK4 , RK45 , RK7 8 ) that are contained in GOPP 
are available to the multirevolution algorithm. STEPR represents a 
new concept for orbit prediction' programs; any orbit predictor can 
have the multirevolution technique as an additional option, thereby 
making it several times more efficient for long predictions.’ 

STEPR has been used as a prototype program and tool for the 
following investigations: 


Notes : 1 . 


2 , 


3 . 


4 . 


5 . 



a. Determine the best set of total energy elements and 
numerical integration methods for applying multirevolution to a 
variety of orbit types (ref. 33). 

b. Determine the accuracy limitations of the multirevolution 
method (ref. 33).- 

c. Study the orbital mo.tion of the proposed solar power 
satellite- over a time span o'f 30 years. 

d. Lifetime studies of near Earth satellite. 

e. Evaluation of refinements to the raultir evoluti on- algo- 
rithm. 

1.7 Concl'udi ng ’ R emarks 

A discussion has been given on several new methods for orbit 
prediction. Included' are numerical, analytica-1 and seminumerical 
methods. Comparison test results have been summarized. Suggestions 
have been made on the use of the new methods for production applica- 
tions. 

For numerical orbit prediction methods, it was found that the 
’•total energy" element formulations produced more efficient results 
when compared to classical formulations. A comparison between the 
two total energy formulations, KS- and PSI*', showed no major differ- 
ence for circular orbits. However, for eccentric orbits, differences 
become apparent because the formulations have -different independent 
variables. It was found 'that, if the independent variable was suited 
for a particular orbit, then more efficient results could be realized. 
Several numerical integration methods were investigated in conjunction 
with’KS and PS<{' formulations'. For circular orbits, all integrators, 
including the fourth order Runge-K'utta (RK4>, proved to be, competitive. 
In eccentric orbits, it was seen that a fixed step Adams method (AD9) 
was the most efficient integr.a.tor , if the formulation's indepen- 
dent variable was suited for -the orbit.. Otherwise, a fifth order 
Runge-Kutta (RK45)' with its excellent step size control proved 
more efficient. 

Accurate and concise analytical solutions have -been obtained 
through' the use of the canonical total energy elements. The solu- 
tions require .a negligible amount ,of computer time and are feasible 
for orbit predictions if the unmodeled forces are small. 

For 'long term, integration , it wa-s found that a seminumerical 
orbit prediction method, based on the multirevolution technique, 
results in a large .savings in computation time. Yet the errors 
in the solution are still smaller than the predicted errors due tc 
force model inaccuraci-es . This versatile method may -be applied to 
any numerical orbit prediction technique. 
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2.0 MATHEMATICAL MODELS OF . THE PERTURBING FORGES' 


2.1 Introduotion 

The environment of an Earth satellite includes a variety of 
forces and effects that actively perturb its orbit from the ideal- 
ized two-body motion. This part of the report will investigate 
these forces and their impact on the accuracy and cost of orbit pre- 
dicti ons . 

A computer program for numerical orbit prediction usually has 
the four basic modules: 

a. Input/output and initialization 

b. Formulation of the satellite differential equations 

c. Numerical integration routine 

d. Mathematical model of the perturbing forces. 


Under consideration here are the various models that are available 
to make -up the fourth module. 

Analytical orbit prediction methods are somewhat diffferent 
from the numerical methods in that they usually do not exist as a 
combina-tion of modules. The formulation, Integration routine and 
force model are contained in the same set of equations. If neces- 
sary, they can be separately. -defined . However, it is not the inten- 
tion to discuss the solution me.thod in this section. The important 
thing is that each orbit prediction method (nu.merical or analytical) 
consists, in some way-, of the above modules. In particular, the 
force model, f once it is defined, plays a large part toward determin- 
ing the applicability of the resulting computer program.- 


■Accuracy is of fundamental importance -in any computer program. 
Sources of possible error in predicting satellite trajectories are: 

a. Roundoff and truncation errors, 

b. Inaccurate initial position and velocity, and 

c. Inaccurate models of the perturbing accelerations. 

Roundoff errors result from the fixed word length that exists 
in all computing ma'chines. That is, -only a limited number of digits 
are carried throughout the> arithmetic operations. Algorithms that 
require more arithmetic operations will be' more affected by roundoff 
errors. From this point of .view, the fewer numerical integration 
steps needed to compute a trajectory, the better. Also, roundoff 
errors' are very machine dependent because of. different word sizes 
and methods of roundi‘ng. 


Truncation 'errors are 
that is carried out on a c 
"discretized” and computed 
and independent variables. 


due to the finite difference calculus 
omputer. Continuous functions are 
at finite intervals of the dependent 
This is mathematically equivalent to 
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truncating the infin 
"truncation error." 
order formulas have 
and less truncation 
The accuracy is usua 
at each step, in ord 


ite Taylor series expression, hence the term 
With numerical integration methods, the higher 
agreement to more terms in the Taylor series 
error. Thus, they are generally more accurate, 
lly controlled to be within prespecified bounds 
er to limit accumulation of truncation errors. 


Section 1.0 discussed truncation and roundoff errors that re- 
suit from a variety of methods; i.e., the numerical accuracy in 
solving the programed equations. 


The second error source, initial position and velocity error 
is related to physical accuracy because, typicaliy, they are ob- 
tained from an orbit determination program that has limited accur 
(due to roundoff and truncation errors, tracking errors, etc.). 
When input into an orbit prediction program, their slight deviati 
from the exact values can cause an appreciable deviation between 
real and computed trajectories. These initialization errors are 
studied in detail in this report. If expected values of the init 
errors are known, their importance relative to roundoff, truncati 
and force model errors can be determined. 


S., 
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This section is concerned with the third error source, inaccu - 
rate and/or incomplete mo_del_Qf. perturbing forces . Generally, the 
models differ in their accuracy and complexity. They must all 
ultimately be based on empirical mathematical formulas and observa- 
tional data. The physical constants (such as mean radius of the 
Earth, gravitational constant, air drag coefficient, etc.) are deter 
mined from such data and are necessarily limited in accuracy, due to 
observational inaccuracies and roundoff in computations. It is not 
the intention of this study to verify or refine the agreement to 
observations. Instead, several existing models are systematically 
compared to determine their efficiency and accuracy when used in an 
orbit prediction program. Efficiency is determined by a model's 
computer execution time and storage. It is desired to find -the 
least complex force model that delivers a required accuracy in a com 
puted trajectory. 

A necessary condition for using numerical orbit prediction 
programs to study force model effects is that the numerical errors 
be less than the force model errors to be studied. The analysis 
and comparisons carried out in section 1 .0 prove that this condi- 
tion is satisfied in the case of the programs (KSFAST and GOPP) 
that are used in this study. 

Typical forces affecting the motion of an Earth satellite- 

are; 


a. Atmospheric drag 

b. Nonsphericity of the Earth 
0 . Sun and Moon gravity 

d. Reference coordinate system inaccuracies 

e. Solar radiation pressure 
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f. Uncoupled attitude maneuvers 

g. Vehicle venting 

An indepth discussion of the first four are given i-n this report. 
Density models of the upper atmosphere are discussed in section 2.2. 
The USSR-ASTP and Jaochia- models are compared. A new dynamical 
(time dependent) density model is presented and tested against the 
Jacchia model. In section 2.3, an analytical Sun-Moon ephemeris 
model is described and compared to the stored ephemeris of the JPL 
tape. The geopotential model is discussed in seeti.on 2.4. It is 
shown how the different terms in the geopotential expression affect 
the trajectory. Suggestions are made on a simplified geopotential 
model. Section 2.5 concerns the effects of noniner-tial coordinate 
systems on the trajectory. Finally, recommendations are given in 
section 2.6 on the correct force models to be used when predicting 
the following types of trajectories: 

Shuttle- tyjie near Earth circular orbits 

Elliptical transfer orbits (e S' 0.7) 

Geosynchronous orbits 

Near Earth orbit lifetime problems 

The complexity of the various forces, and the diversity of 
satellite shapes and orbits, makes it very difficult to consider all 
possible perturbations. This study has,, of necessity, been somewhat 
limited in scope. The last three of the above mentioned affects will 
not be discussed. However, it will be made plear in each subsection 
as to the scope and limitations of each force model studied. 

2.2 Atmospheric Density Models 

Atmospheric drag is an important perturbation on the orbits of 
near Earth satellites, such as the shuttle orbiter. The drag force 
model includes two effects: 

a. The aerodynamic interaction of the satellite with the 
upper atmosphere, and 

b. The variable atmospheric density. 

The first depends on vehicle characteristics; i.e., the satellite's 
shape and body attitude. The second depends on altitude, time of 
day, geographical position and date- 

The acceleration a due to drag is based on the. empirical 
equation 

( 2 - 1 ) 
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where p is the atmospheric density at the position of the 

satellite, 'V 'is the velocity magnitude and direction of the 
satelli-te relative to the atmosphere., and is the drag coeffi- 

cient. The ballistic number B is defined as the weight divided 
by the cross sectional area of the satellite. 

It has been determined experimentally that the drag . coefficient 
of a vehicle moving in a rarlfied atmosphere is approximately 2.2. 
Therefore, the studies described in this section ha.ve been based on 

Cd = 2.2 

The ballistic number depends, in general, on the angle .between 

-f. + •+ 

the vectors V and h, where si is the body orientation vector., 

For the shuttle orbiter, B can vary between (approximately) 50 and 

400 lbs/ft2. For this study, an intermediate value of 

B = 100 pounds per square foot 

has been used throughout; i.e., B was held constant in all com- 
parisons . 

2-.2.1 Description of the models. .- The density at any point 
above the Earth’s surface generally depends on altitude, time of day, 
and level of solar activity. These variations, if incorrectly model- 
ed, may result in inaccurate orbit predictions for satellites that 
are near the Earth. Nondynamic models such as the 1962 U . S . Standard 
Atmosphere do not include the time depend ence and solar activity 
effects. Four density models were investigated (refs. 34, 35, 36, and 
37)*. 

a. An analytical nondyhamical model based on the 1962 U.S. Stan- 
dard Atmosphere. This model is described in reference 38 and 

will be called the AMDS model. It uses exponential functions to 
describe the 1962 atmosphere. In a trajectory prediction pro- 
gram, it uses a. negligible amount of computer storage and execu- 
tion time. 

b. The Jacchia model (ref. 39). This model is composed- of 
two parts: 

(1) The determination of the exospheric temperature as 
a function' of position, time, and the solar and geomagnetic ac- 
tivity. 

(2) The determination of density as a function- of the 
exospheric temperature and altitude. 
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In addition to the solar and geomagnetic activity, the model con- 
tains semiannual and diurnal atmospheric variations. Empirical 
formulas based on these variations are used to obtain the exo- 
spheric temperature. 

c. The USSR (Russian) model used in the Apollo-Soyuz Test 
Project (ASTP) mission. It includes the four effects mentioned in 
(b) above and is based on the tracking data of Cosmos satellites 
for the period from 196^1 to 1970 inclusive, , It is not as versatile 
as Jacchia's model since the solar radiation intensity data is 
assumed fixed for the year 1975 (ref. MO). 

d. A new analytical model (AMDB*) which contains diurnal 
variations (refs. 36 and 37). In addition, it can be initialized 
to agree with Jacchia's model. It requires less computer storage 
and runtime than Jacohia or USSR. 

The required input data for each model is_ as follows: 

AMDB - satellite altitude above the Earth 

Jacchi a - satellite altitude above the Earth, 

- position vector of the satellite relative to the 
Earth, 

- position vector of the Sun relative to the Earth, 

- Julian date of the, satellite ' s ‘ state vector 
Russian - satellite altitude above the Earth, 

- posi-tion vector of the -satellite relative to the 
Earth , 

- position vector of the Sun relative to the Earth, 

- number of hours from midnight December 31 
AMDB * - satellite altitude above the Earth 

- position vector of the satellite relative to the 

Earth , 

- Julian date of the satellite's state vector 

The AMDB* model was developed by Gus Babb (NASA/JSC), 

Stephen Starke (ACM) and Alan Mueller (ACM) when it became appar- 
ent that the three other models had significant deficiencies. The 
most notable were; 
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AMD.B - The exponential formulas were not accurate en'ough. 

In some cases,, this model gave no better results 
to orbit predictions than' if air drag had been 
completely neglected. 

Jacchia - This model involved large costs in terms of computer 
storage and execution time. 

USSR ■ - Poor results were ob'tadned' for epochs other than 

1975 . Also, it was grossly inaccurate at 'the lower 
altitudes (120 km). 

2.2.2 Method of evaluation .- The purposes' of thi-s study are . 
twofold: (1) determine, the effects of atmospheric density varia- 

tions on the path of a near Earth satellite, and (2) evaluate the 
various density models in terms of accuracy, and efficiency. 


It would be' desirable to compare each model directly to the 


actual density of the upper atmosph 
This, is not possible because direct 
available. In fact, the density is 
effect on satellite motions. This 
the Jacchia and USSR models. There 
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Since both the Jacchia and USSR models' were developed indepen- 
dently, arid' .since both are based on sate'l-lite tracking data, a con- 
servative estimate -of thei-r accuracy -can be determined by comparing 
them to each other in epoch 1975. Output densities are compared in 
table II. 


•TABLE II.- DENSITY OUTPUT, JACCHIA VERSUS USSR 


Altitude , 



120 

240 

320 

50.0 


Di f f erenc e , 
percent 

71-.2 

8.9 

1 .4 

9.8 


The large difference' at 120 km is due to the inaccurate USSR- 
model at lower altitudes.. It therefore should not be used when the 
satellite's '‘altitude is less than approximately 200 km. 
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Trajectory, comparisons .were made' by -exami-ning' the predicted 
posi-ti'on of a* near Earth satellite '(‘hg = 220. 1cm, hp = 380 km). 

A complete geopoten.tial model wasi used.. Predic.ted position based 
on Jacchia and USSR density models are compared in Table III for 
epoch- 1975 . 


TABLE III.- POSITON DIRF^ERENCE, J-ACCHXA V.ERSUS. USSR • 


Prediction 

Position • 

interval, da-vs 

difference, km 

0.5 

0.,4 

1.0 

• 1-.7 

5.^0 

‘3-9 . 0 


The Jacchia model is considered to be more complete than the 
USSR model, because i-t can be applied .to the full- range of altitudes 
and epochs (ref., 34). Aiso, it is ■ in . wide -us^e .-at NASA/JSC. It was 
used,, therefore, as the -reference for these studies. 

2.'2.3 Effects of the_Sun on the upper, atmosphere .- In addition 
to the nearly exponential variation with altitude, the upper atmo- 
spheric density varies according to the time of day (diurnal), time 
of year ( seasona.l ) , -.and I'evel of solar activity.. The .diurnal effect 
is a “bulge" on the atmosphere caused by solar heating of the sunlit 
side. A change in ^seasons causes the bulge to change latitudes. 

It has been found that these three effects can cause a signifi- 
cant variation in the density at any point of altitude 500 000 feet 
or higher. This section shows the variation as a function of altitude 
of the diurnal and, solar activity effects. 

The J.acchia model was used to determine the atmospheric 
density at altitudes between 120 km and 600 km. Figure 3(a) 
shows how the density varies between these altitudes. The four 
curves show the difference between the atmosphere in sunlight 
and darkness, as well as the effect of solar Intensity in 1970 
and 19-75. The Sun shadow curves show the diurnal effect. The 

reference density is Oq = 1.225 kg/m^ . 

! The curves in- figure 3,(a) converge at about 500.000 feet 
(15'2 km). A.t. this al.ti.tude and below .there-- is. good, '.agreement 
between Jacchia and the AMDB exponential model. 
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More explicit details of the diurnal variations are shown 
in figure 3(h) where the density ratio function -in p/Pq is 

plotted against the Sun's hour angle. The daily variation is 
apparent . 

To determine the effects of these variations, satellite orbit 
predictions were carried out. using the Jaochia and AMDS models, and 
also for the case where air drag is neglected. The orbi-t, conditions 
were the same as for the case in section .2.2.2. Positions in the 
orbi't are- compared against, the Jaechia-based solution i'n table IV . 

TABLE IV.- POSITION DIFFERENCE., J'ACCHIA VERSUS AMpB,. NODRAG 


Prediction 
interval, davs 

Position 

difference AMDB . km 

, ■ .Position 
difference NODRAG 

0.5 

11.1 

6.7 

1 .0 

l|4 .8 

27.8 

.5-.-0 

1127.0 

,64 9.0 


For this particular ease, the nondynamic exponential, model 
gives no better results than would' be obtained from completely 
neglecting atmospheric drag. These results lead to the conclusion 
that any accurate orbit prediction program must make use of a den- 
sity model that includes the three time dependent effects. 

2.2.H A new analytl_cal. at-mjj.s pheric density model.- Inspection 
of curves in figure 3(a) suggest that, at any given time, the atmo- 
spheric density function (ordinate) may be represented by 

-An(P/Po) = F(-z) (2-2) 

where 

P = atmospheric density j^g 

p^j = reference density at sea le.vel (1.225 — =) 

. ■ 

z = altitude above the oblate Earth in kilometers, 

and F(z) is a rational polynomial. The expression for p 
would be ' 
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(a) Variations with altitude and epoch . 
Figure 3.- Atmospheric density variations. 



(b> Variations with position of the Sun at 500 km altitude 
Figure 3.- Concluded, 


The coefficients in F could be determined to give agreement to 
the curve that corresponds to a given epoch. This procedure would 
account for seasonal and solar radiation intensity variations,, at 
least over a limited time of a few days or weeks. In addition, it 
was thought' time dependent trigonometric terms might be added to 
F in order to include the diurnal variations. 

This approach has -been used to derive a new analytical at- 
mospheric den sity . model (refS'. 36 'and 37). F in equation (2-3) 
is represented by 


F = a,,+a„z+^+B. 
X z z 


(2-J{ ) 


where 


B = -b (z-152) eos (>|)-35°) cos 8 (2-5')-- 

, ag, a^. and b are constants to be determined. The function 

B effectively simul.ate.s the diurnal bulge, of the. atmosphere. >l> 
and 8 gi-ve the angular distance of the subsatellite_ point from 
the center of the bulge. Given the right ascension and deolinati-on 
of the Suii (ogj-Sg) and the, right ascension and declination of the g. 
vehicle (ay,6y)i then 


8 - - «g cos 4i 

The constant of 35° in equatio.n ( 2-5 ), represents .the, fact that the 

bulge lags 35° behind .the Sun, as viewed from the surface of the 
Earth. 

The vai'ues of a.| , ag, a^, are calibrated for agreement - of ■ 

p with the va-lue predicted ,by the Jacchia model'. Three- points are 
chosen; z-| = 152 km, Z 2 = ^00 km, = 600- km. They are positioned 

over the surface of the Earth such that B = 0-, i.e. i]; - 35 = 

B = I . f.,.} f^, ^3 SI’S the -values of 



as determined from the Jacchia model. Then a.| , a 2 , a^ are com- 
puted sequentially from 
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^ ^2 h ^ 1^2 


(2-7) 


®3. 

^1 “ ^1 Vl " z^ 


To compute b, choose z = 400 1cm and ' ^-35° = I80®j ‘ g = O' 
f is obtained from- the Jacchia, such that 


b-(£-«l- 


^2^ ~ z 



(z - 152) 


( 2 - 8 ) 


Equation (2-5), and hence the AMDS* model-, is valid for 
z 2152 km. For altitudes b'elow T52 km, the AMDB model should be 
used since it is in close agreement with Jacch-ia for that a-ltitude 
range (see fig. 3(a)). 

It was expected that this model would give good agreement 
with 'the' Jacchia moUel ‘but require less ’c-omputer storage a-nd exe- 
cution time. To inyestigat.e this concept, several orbit prediction 
experiments were carried out. 

2 , 2.5 Orbit predictipn experiments- .- Numerical orbit predic- 
tions were carried out using the GOPP program (ref, 25), The orbits 
used in the comparisons are shown in table V, 


TABLE V.- ORBITS USED IN PREDICTION EXPERIMENTS 
■ ^ 



A 

B 

C 

D 

E 

F 

perigee ('km) 

220 

300 

300 

220 

Same as 

166 

apogee (km) 

380 

600 

600 

380 

A except 

433 

eccentricity 

.012 

.022 

.022 

.012 

epoch 

.02 

period (min) 

90.5 

93.6 

93.6 

90.5 

is 

90.5 

argument of perigee 

0 

0 

180° 

0 

1.2 : 00 

0 

ascending node 

0 

0 

0 

0 

January 1 , 

0 

inclination 

epoch' 

30 ° 

12:00 

3qO 3qO 

January. 1 , 

90° 

1975 

1977 

30 ° 

Same as 


A ,B , C , D 
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Results of the orbit prediction experiments with different density 
models are shown in tables VI(a) through 'VI{f). All comparisons are 
with respect to the Jacchi’a model- The "No Drag" case_ is also inc'luded- 


TABLE VI.- POSITION DEPENDENCE ON DENSITY MODEL 

(a) Orbit A 


Time of 
integration , 
davs 

Ppaltipn 

difference, km 


No dr.aR 

RUSSIAN 

AMDB» 

0.5 

8.9 

0.2 

0.3 

1 .0 

36. 1 

0,8 

1.5 

5.0 

895 

16 

42 



(b) 

Orbit B 

" 


Time of' 
integration , 
davs 

- 

Position 

difference, km 


No draK 


RUSSIAN 

AMDB* 

0.5 

.83 


.03 

.07 

1-0 

3.16 


-.07 

.28 

5.0 

81 . 5 


^1-0 

6.0 
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TABLE VI 


Contlriued 


(e) Orbit G 


Time of 
integrati 9 n , 
da vs 

No draE 

Position difference, km 
I^U;?SIAN 

AMDB» 

0.5 

.5 

0 

.01 

1 .0 

1 .8 

.02 

.04 

5.0 

'*7. 5 

.8 

.5 


(d) Orbit D 

Time of 




integration , 

Position 

difference, km 



No draE 

RUSSIAN 

AMDB* 

'0.5 

8.6 - 

. 1 

.3 

1 .0 

34.4 

.5 

1 .5 

5.0 

836 

10.3 

44.9 



(e) 

Orbit E 



Time o‘f 
integration , 
davs 

No draE 

Position 

difference > km 
RUSSIAN 

AMDB* 

0.5 

13.9 


3.4 

0.4 

-1.0 

56.0 


19.2 

1.9 

5.0 

1404 


494 

■ 41.6 


39 







TABLE VI. - Concluded 


(f) Orbit F 


Time of 
integration , 

P.osition . 

difference, km 

davs 

No draK 

AMDB» 

0.5 

40.9 

.4 

1 .0 

164. ‘4 

2.0 

5.0 

3944 

25.9 


Notes: 1. The position difference is almost entirely in the • 

downrange direction. The largest effect of drag is 
to perturb the position, of the satellite in its- orbit. 

2. All of these cases show that drag can cause a signifi-". 
oan.t perturbation for orbit predictions of one day or 
longer. 

3. The new analytical model CAMDB*) shows goo'd a'greem'ent 
with^ Jaochia, although not quite as good as USSR for 
epoch in 1975. 

4. For epochs other than 1975 (orbit E) the USSR model is 
not very accurate. This shows the effect of the solar 
activity level. The USSR. model is based on. the solar 
activity in 1975. 

5. AMDS* has about the same accuracy, regardless of epoch.,' 
(Compare orbit A with orbi.t E.) 

2.2,6 Addi'tional remarks .- It has been shown that atmospheric 
drag can cause a" large perturbation of a satellite's position in 
its orbit. In addition, time dependent variations in the atmo- 
spheric density ca"n also cause large perturbations. Therefore, 
orbit prediction programs for near Earth orbits must include a 
dynamic model of the atmosphere. 

The Jacchia density model has been compared with the USSR 
model. They show good agreement for certain orbits. The USSR 
model, however, is not valid for altitudes lower "than 200 km or 

epochs other than 1975^. Also, both models are rather ineffi- 
cient in terms of computer storage- and runtime. ' ' 


^If given the correct solar activity data, the USSR model 
could, perhaps, be calibrated to work for any epoch. 



A n.ew analytical model (AMp’B*;) has b'een presented. This 
model 'Was developed jointly by, Gus Babb ( NASA /-JSC'). > Stephen Starke 
(ACM) and Alan Mueller (ACM)* ’ It can be calibrated to agree with • 
the Jacchia modeL at any given .epoch.' The chief advantages of this 
model are its small computer storage and 'execution 'time requirements. 


Execution 
are shown 
gramed -on 


tlme^ and storage requirements -for the four models studied 
in table VII. This data is based on the models being pro- 
the Univac 1110 system. 


table VII.- EXECUTION TIME AND STORAGE REQUI-REMENTS 
(.Atmospheric Density .Models)- 


Model 

Jacchia 

USSR 

AMDB 

AMDB* 


Eveoution time, ms 
1.6 
2.6 
0.2 
'0.3 


Storage, words 
1252 
'57-5 
<100 
3=11 . 


It can be' seen that AMDB*' offers considerable advantages, in 
speed and storage over Jacchia and USSR, and is almost as fast as 
the simple exponential model (AMDB"). This increased speed is most 
useful for long orbit predictions, lifetime studies, or in appli- 
cations where execution time must be limited, such as in an onboard 
computer . 


2.3 Geopotential Model 

The mathematical model for' the gravitational potential of the 
Earth is given as a function of the satellite’s position with re^ 
spect to an Earth fixed coordinate, system. This function is the well 
known solution of Laplace’s equation in terms of spherical harmonics. 
The perturbing accelerations are obtained as the partial derivatives 
of the geopotential. They can be computed using nonsingu-lar recursive 
equations (ref. ) . 


^Some additional savings in execution time for the USSR could 
be realized if single precrsi-on programing were- used. 
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The figure of the Earth is described by the numerical coef- 
ficients in the geopotential expansion'. These coefficients have 
been computed from satellite observations and- surface gravity 
measurements (ref. 42). Although they may sometimes be orbit de- - 
pendent, the thorough analysis of' reference 42 'provides coeffi- 
cients that describe the true Earth to a high precision. These 
coefficients were taken as the reference model for the comparisons 
described in this section. 

Orbit predictions based on the fully determined geopotential 
will give the best accuracy but can be prohibitively expensive in 
terms of computer execution time. ,It is des.ired, therefore, to 
determine the effects of neglecting many of the higher order terms. 
Attempts Mere made in references 43 and 44 to determine satellite 
position errors resulting from predictions based on a truncated 
geopotential . The numerical results showed that any neglected term 
could produce a linear error growth. However, the analytical sat- 
ellite theories and observations (table 1, ref. 42) indicate that 
the tesseral terms produce only periodic perturbations of small 
amplitude (200 nautical mile altitude). Additiojial numerical ex- 
periments described in this section showed that, when properly ini- 
tialized, the linear error growth (due to neglected tesseral terms) 
in the predicted orbit could be avoided. 

2 , 3.1 Numerical pyperiments . - Studies were conducted comparing 
a Jg perturbed solution with a Jg *^22 perturbed solution. 

CJ 2 being the second zonal harmonic and *^ 22 ' the second tesseral 

harmonic.) The J 2 perturbing potential is a function of only the 

position magnitude and latitude, while the J 22 potential is a 

function of the .position magnitude,, latitude, and longitude with re- 
spect to the Earth. Two different cases were run for the J 2 versus 

J 2 + J 22 comparison. All were given the same inertial initial 

conditions. However the initial longitudes with respect to the 
Earth were different; i.e., different Initial epochs. The initial 
inertial conditions were: 


6 6.7 7 . 7 6 6 km . 

Q = 

0 

0 

0 

.015 

(j = 

0 . 0 ° 

30° 

M = 

0 . 0 ° 


In case 1 the hour angle was = 104.8° and in case 2, 

The position differences between the J 2 3*^^ *^22* 

shown in table Vlll(a) for both oases as a function of time. 


Gq = 149.9°. 
solutions are 



As one can see from table' V-III(a), 'case .1 appears to show a 
secular position difference trend stemming from the J 22 

turbing force. While in case' 2 th'e 'position difference seems to 
be periodic and much smaller. Case 2 verifies what has been shown 
analytically and from observations.. 

Another numerical study was made comparing an- eighth order 
zonal model with an eighth order eighth degree zonal and tesseral 
model. Again two different cases for the hour angle were chosen. 
In ease 1 the hour angle = 42.9 while in case 2 s-q = 100.6. 

The initial conditions were the same as those in reference 44. 


6629.656567 km. 

fl = 

4'5° 

. 01 

w = 

45° 

45° 

M = 

45° 


The Julian date of 2442332.17 correspnds to the hour angle in 
case 1. The fe'sults are. shown in table ‘Vlll(b) . . 

The resui’t's in case 1 , table VTII'(b-) are similar to- those 
shown in reference '44'-. Again, table Vlll(b) shows- that the error 
growth is strongly dependent on the initial epoch-. 

TABLE VIII.- POSITION DEPENDENCE OK GEOPOTENTIAL MODEL 

(a) 2x0 versus 2x2 



Position 

difference.' km 

dfl ys 

Case 1 ' 

Case 2 

0.25 

1 . 35 

0.43 

.50 

4.91 

.005 

1 .0 

4,07 

. 02 

1.5 

6.09 

.03 

0 

CM 

19.54 

.04 


43 


TABLE VIII..- Concluded 


Time 

- days, 

.25 

.50 

1.0 

1.5 

2.0 

2.3.2 Theorei.loal background .- The basic reason for the 
results of the two previous examples is that the right-hand side 
of a differential equation cannot be arbitrarily truncated without 
possible drastic effects on the solution. Consider the .single or- 
dinary differential equation for x, 

X = F(x,t) (2-9) 

A similar differential equation for y is 

V = F(y>t) + eG(y,t) (2-10) 

where le|<<1 and G is periodic in y and t. Let x*(t) 

and y*(t) be solutions to equations '( 2-9 ) and (2—10), respectively 

such that 

x*(t^) = y*(to) 

Define 

^(t) = x(t) - y(t) 

Then 

5 = F[c + y(t),t] - F[x(t) - - eG[x(t) - C«t] 

When ? is initially zero, 

? = F[? + y«(t),t] - P[x«(t) - C,t] - eGtx»(t) - C,t3 (2-11) 

where x*(t) and y*(t) are known functions. According to the 
theory of differential equations, there is no reason to expect that 


4M 


(b,)- 8x0 .versus' 8x8 



2 : 45 . 

1 .29 

7.3 

.31 

16.49 

' .92 

24.72 

1.84 

28.9 

2.02 



the solution. to e.quation will remain sDia'll.. ■ There "is thus 

no theoretical justification, for arbitrarily neglecting geopqtential 
terms . 

Much is known about the solutions of the differential equa- 
tions in celestial mechanics. As mentioned earlier, the tessera-1' 
terms produce periodic effects in the solution'.. Therefore, an 
average over one revolution can be used to determine a "mean" 
mean motion that will compensate for the neglected tesseral term.. 

The mean mean motions in case 2- (tables Vlll(a) and 7(b)) agreed, 
more closely than in case 1 . 

2 . 3.3 Comments and recommendations .- Based . on. -the results 
of these studies and those documented in references and 44, the 
following recommendations can be made: 

• For very accurate predictions of satellite, position (near 
Earth orbits) the full geopotentiai model ,inust be used. 

• Predictions of orbit size and shap.e can be based on a 
zonal model only. 

•, Mean initial elements must be used before any tesseral 
term can be neglected.. 

•Practical realiza.tion of the third’ recommendation requir.es 
additional work in the following areas: 

• Development of a rigorous theoretical justification for 
using initial mean elements with' a simplified geopoten.tial . 
This .should, .make, the- connections between analytical theories 
and numerical solutions. 

t Development of a fast numerical routine -for computing ini- 
tial mean values for the elements. Impact bn existing pro- 
- -grams- should .be minimized with th.e, goal of minimum storage 
and. minimum execution time for the ini.tial"ization . 

• . Development of an efficient m.e.thod for .inclu.ding important 

resonance effects. It is known that the motion of the 
satellite can be significantly perturbed by, certain res- 
lonant .geopotenti-al terms. These canno.t be neglected .j. 



me gravitational forces of attraction due to the Sun- and the 
Moon can have an important effect on the orbit of a satellite, par- 
ticularly. for hi-g-h altitude orbits. Given the mass and position of 
a perturbing body, one may calculate its perturbation on the orbit-. 
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Let X be the position vector of the satellite, referenced to 
some inertial coordinate system centered -at- the Earth." Then its 
acceleration is 


-y u -> ^ • 

x+-x = T;, +F_+R 

r 1 © 


( 2 - 12 ) 


where 
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X - r^ , 
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and F 

0 

are the 

disti 


(2-13) 


■© 


Sun, respectively, and E represents any addTtion'a-1 accelerations. 

r and r are the positions of the Moon and Sun, respectively, 

» 0 ■ . ■ _ • s 

referenced to the- Earth centered inertia-1 coordinate system. The 
"mass ratios of the Moon and Sun are known -to a~-very "high precision; 
The gravitational parameter ’is 


,U = M, 


where k is the gravitational constant and 
the Earth ' 

The first term in. 'the brackets of- and 

"direct" term and represents the gravlta-tional a-ttraction of the 
satellite by the' perturbing body. The second term-, called .-the "in- 
direct" term, represents the perturbation of the motion of the 
Earth. Equation (2-13) assumes point masses for the Moon and Sun. 

+ 

R in equation (2-12) represents perturbations such as atmospheric 
drag, nonsphericitjr of the Earth, "etc. 

A'ccording to equation (2-13), the ‘luni-s'olar perturbing ac- 

celerations ^ and F can be directly computed, provided that 


-is- the ‘mass of 

is called the 

CO 


46 



numeri’oal integration of equation 

the position of the Sun and Moon 
Two methods for obtaining this 
d and are discussed in this sec- 


a. The stored ephemeris data on the Jet Propulsion Labora- 
tory (JPL) tape. A polynomial interpolation is used to obtain 

intermediate values of and Tq. This ephemeris is considered 

to. -be- the most accurate .availahle, and will be used as the reference 
for this study. It was derived based on radar observations and sat- 
ellite tracking data.’. . . 

‘b. Analytica-1 formu-las that give *the elements of the Sun and 
Moon as a function of time. P_ositi-ons can -then be computed from-, 
these elements-. This type' of ephemeris is sometimes preferred' be- 
cause it does not require the tapei read and. table lookrup operations, 
However, the analytical formulas are ‘generally not as'’acourate as 
the JPL tabulated data. 

The lunl-solar perturbatons will' be shown for 'near Earth, geo- 
synchronous and elliptical orbits. An analytical ephemeris model 
will be discussed and compared to the JPL model. Of particular in- 
terest are their accuracy and ' their computer storage and runtime 
requirements. Finally, suggestions are made for better utilization 
of a stored ephemeris. 

'The Initial parameters of .the three orbits studied are shown 
in -table IX . 


TABLE IX.- OBBITS USED IN PREDICTION EXPERIMEN-TS 

(Luni-Solar ) 


X, ^ and r^ are known. The 

(2-12) requires, therefore, that 
be known as a function of .time, 
information have been investigate 
tion: ■ 


.Geosynchronous, km 

Near-Earth", km 

Eccentric, km 

35 864 

.220 

1 000 

35 864 

380 

39 462 

0.0 

0.012 

0.723 


I 


.0 


30 


30 


TABLE IX.- Concluded 
(Luni-Solar) 


Period 


geosynchronous, km , Kpi E ggg ntr?-g> km 


2H hrs 


90.5 min 12 hrs 


E poch 


Noon January 1, 1975 


2.4.1 Magnitude of luni-solar,- perturbations ■ - The total 
effects of the Sun and Moon on the three orbits was determined by 
comparing the output position vectors, of - two cases of orbit pre- 
diction. In one case, the luni-solar gravity per turbati ons. were 
included using the JPL stored ephemeris (ref; 45). Luni-solar 
gravity was set to zero in the other case . Results of the com- 
parisons are shown in table X(a). 

The individual effects of Sun and Moon gravity are shown in 
tables X(b) and X(c), respectively. 

TABLE X.^ LUNI-SOLAR GRAVITY EFFECTS ON POSITION 
(a) Combined Sun and Moon 


Prediction 
interval, davs 


Near-Earth, ra 


0.5 

3.7 

3 6;. 8 

27.3 

1 

11.1 

75.7 

54.4 

5 

57.5 

39K5 

174.7 

30 

281.0 
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. TABLE .X,..- Conc'luded* 
Cb). ’Sun’ on'.l'y 


Prediction 
interval, da-vs* 

0,5 

fienRvnohrnnous km 

:Nea r.-Earbh . m , 

Eccentric, km 

2.3 

= - '31- -3 . 

- 2.5 

1 

' 5.0 

62.9' 

4'. 7 

5 

25-. 0 

321.5 

1-4.9 

■ 30 ' 

' 152.4 

- 



(c-) Moon 

only 

• 

Prediction- ' 
interval, days 

0.5 

Geosynchronous, km 

'Near-Earth, m 

Eccentric, km 

4.2 

5.8 

•24,. ,9 

1 

6.1 

13.1 

49.8' 

5 

32 .. 7 

70.8 

. .174 .-8 

30 

128.7 

— 



Notes: 1. The errors in each case have a linear growth. This indi- 

cates -that the decision to. include lun’i.-solar -per turba- 
tions in an orbit prediction depends on the- expected pre- 
di'c.t-ion inbe'rva'l a'hd requl-red -accuracy. 


2. The Moon has a much -larger e'ff.ect. than- the Sun on the 
eccentric orbit. On the other handf the Sun has a 
larger effect on the near Earth- orbit. 

3 . Compared to air -density uncertainty effee.t's (table IlDi 
the- luni-solar perturbations on a near Earth satellite 
are smalT. 


The linear growth irid-i-cates- that the source of the error 
iS' a slightly wrong value of the "mean" mean motion, simi- 
lar to t'hat' discussed for the g‘eopotential in section 2,3. 
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2.4.2 Analytical mod el versus JP L_aj;_&red data .- An alternative 
to the JPL stored ephemeris table has been investigated in references 
46 and 47. It is based on expressions that give the mean orbital 
elements of the Sun and Moon as functions of time, and it is therefore 
referred to as an analytical model. 

Expressions for the Sun were obtained from the American Ephermeris 
and Mautical Almanac (ref. 48) and are documented in reference 46. 

Brown's lunar theory was used for the Moon's analytical ephemeris (source 
was ref. 49) and the mathematical expressions are given in reference 47. 
The algorithm for these analytical models will be referred to here as 
"ANAL." The JPL data (and associated algorithm for interpolation) will 
be referred to as "JPL". 

The JPL AND ANAL output positions of the Sun and Moon are compared 
in figures 4 and 5. The differences are relatively small and periodic. 

However,, it must be determined how these differences affect an 
orbit prediction. This was done in references 46 and 47, and the 
results are summarized in table XI for the three types of orbits 
shown in table IX. 


TABLE XI.- POSITION DIFFERENCE ,' JPL VERSUS ANAL 



Notes: 1. It is assumed that JPL is the more accurate ephemeris 

model since it is based on precise optical observations, 
spacecraft tracking and radar data. 

2. The errors caused by using ANAL for orbit predictions 
are extremely small and are negligible for many applica- 
ti,ons,. 

3. The choice between JPL and ANAL must be made on the basis 
of operational considerations, such as execution time 
•and storage limita>tions,.- 

Computer r.untime and storage comparisons, are shown in table XII. 

While there is very little difference in execution time, JPL requires. 
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on the average, 4 seconds to initialize and read the data tape. 

In addition, for. Uniyac 1110: demand users., there can be a long de- 
lay while the computer operator finds and loads the magnetic tape. 

The demand mode also requires special procedures for the tape read 
during execution. Computer storage 'is about 50 percent more for JPL. 

TABLE XII.- COMPUTER HUMTIME AND STORAGE, JPL VERSUS ANAL 


Mo-del 
' JPL 
ANAL 


2.4‘,.3 Additional suggestions .- ‘Improvements can be made_ to 
the. JPL algorithm that, may substantially reduce the disadvantages 
mentioned in .the previous section. The suggested approaches are; 

a. Strip the Sun and Moon data from the JPL data set. 

Carryout the interpolation calculations for the Sun and Moon only. 

b . Place the Sun-Moon data 'in a more' easily accessible mass 
storage device, such as Fastran files. This would require seg- 
menting' the data according to epoc'h. 

e. . Use Chebyschev polynomials as the interpolating functions. 
This would provide a bound on errors and require less data storage. 

2,5 Time and Coordinate System Models 

The direction in space of the Earth’s rotational axis is not 
fixed but has two separate motions, called precession and nutation. 
It is usually desirable to express'fhe geopotential accelerations 
in an .Earth fixed coordinate system. These accelerations must then 
be converted to an inertial- coordinate system in order to carryout 
the numerical integration of the satellite differential equations. 

The Earth's rate of rotation is slowly decreasing.' This 
effect, if uncorrected, can introduce errors in time. In order 
to correctly compute the perturbations- due ' to the nonspherical , 
rotating Earth, the Earth's hour angle must be accurately known. 


®This is the computer cost on Univac 1110 for one evaluation 
of the algorithm for Sun and Moon positions. 
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This section will discuss studies that were carried out 
in order to investigate the most efficient methods for including 
the effects of precessioni nutation and Greenwich rotation in 
an orbit prediction program. Also to be discussed are numerical 
experiments, that show the errors that will occur when these 
effects are neglected. 

2.5.1 Di sou salon .- Precession is the steady drifting of the 
Earth's polar axis of rotation around the surface of a cone whose 
axis is the ecliptic pole and whose semicone angle is the obliquity 

(23®27') (see fig. 6). Period of revolution around this cone is 
26 000 years, which corresponds' to a drift of approximately 20 arc- 
seconds per anum at the Earth's pole. Suppose an inertial reference 
system is defined such that the Z-axis i-s in the direction in which 
the Earth's mean pole pointed on January 1, 1950. Then the Earth's 
pole in 1975 will differ from the Z-axis of the basic reference 
coordinate system by approximately 514 arc-seconds. 

Define the precession matrix, P to give the transformation 
from mean equator and equinox coordinates of some reference epoch 
(Basic Reference Inertial Coordinate System) to mean equator and 
equinox coordinates of date. The expression for this matrix 
is given in reference 50. 

Nutation refers to the small os cilia tory ''wobbling of the 
Earth's polar axis around, its mean precesslng position. As a re- 
sult, the value of the obliq.uity of the' ecliptic oscillates about 
a mean value. This transformation represents the difference' be- 
tween the position of the true celestial pole and the mean celestial 
pole. Nutation may be broken up into ,a series of short period term: 
(up to a total of 139 terms). Some of these terms have periods as 
small as 5 1/2 days. The two largest terms correspond to the 
Earth's pole traversing an ellipse of semiaxes 9". 2 and 6". 8 
once in 18.6 years. The six largest of these terms were included 
in the study. 

Thus, the nutation matrix N makes the transformation from 
mean equator and equinox coordinates of date to true equator and 
equinox coordinates of date. .Explicit expressions for, the matrix 
N as a function time are- given in reference 50. 

Greenwich rotation is the daily rotation of the Earth around 
its own spin axis. The rate of rotation of the Earth was derived 
by Newcomb in 1895. Hi's equation became the. definition of "non- 
uniform" time commonly called universal time, UT1. - The amount of 
rotation of the Earth is given in terms of Ufl rather than the 
angle called the Greenwich mean sideral time. in the analysis, 
Greenwich mean sideral time is obtained from UT1. The explicit ex- 
pression for the Greenwich rotation (matrix G) is given in refer- 
ence 50. 
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Figure 6.- Precession nutation, and Greenwich rotation. ■ 



To summarize this discussion, the entire transformation from 


basic 

reference 

iner 

tial coordinates 

to Earth-fixed 

CO 

ordinates, 

and vice versa, 

can 

be given by the 

matrix product 

GNP 

as follows: 


>'EF = ° 

N P 

’"IC 




or 








ric = P^ 

n'I’g-t 

^EF 




where 

r-£^ and 

■"EF 

are vectors in 

the h.efe.rence 

coo 

rdinates and 


Earth-fixed coordinates, respectively. 

2.5.2 Humerical experiments -T'lfe' expressions for matrices 
•G , N, and . P are lengthy, ootit-ainin'g many trigonometric functions 
It is, therefore", desirable when numerically predicting satellite 
orbits to minimize the number of times these matrices need to be 
calou'la'ted'.. A study* was* carried''out to determine the effect of pre 
cession, nutation- and nonuniform .Greenwich rotation on a near Earth 
satellite orbit and to determine the . o'ptimum way to implement these 
effects in orbit prediction, programs.' 


Two 

N, and ,P 
per turbin 
In the se 
at the ep 
is change 
GMP and 
due to pr 
ferred to 


methods of approach were taken.’ In one, the matrices G, 
are computed .-at-each .function evaluation (calculation of 
g' forces).- This method will., be .re ferred to here as GNP . 
cond method of . approach'.,’ these matrices were computed once 
och of Initialization and held’.' constant’- (till epoch time 
d, if ever); This method,-„will’ be referired .to here as GNP 
GNP*’ were • compared with neglec-ti.^ng totally the effects 
ecession nutation and no'nuni form' Greenwich rotation (re- 
as GNP»*); 


In GNP**, the tranaf ormaMon ■ f rom basic reference to- Ear.th- 
fixed system and vice versa’- is accomplished' . by 


EF 


= R 


uni for 




Epoch ^'^IC 

r ”1 'T T 

ric = G N PJ Epoch ^uniform ’"EF 

where ^uniform ® rotation matrix giving a uniform rotation of 

the Earth. 



cos (Ogt 

sin Ugt 

0 

p 

uniform " 

-sin “Et^ 

cos Ugt 

'o 


0 

0 

0 
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is the rate of uniform' rotation of the Earth at initialization 


U£ 

epoch, and [G N is the exact transformation matrix from 

basic reference to Earth-fixed system at initialization epoch. Thus 
u)g and Q = G N P are computed once and ^uniform computed at 

each function evaluation. 

The example considered has the following initial conditions: 


Height of perigee 

Height of apogee 

Eccentricity 

Period 

Inclination 

Epoch 

Perturbations 


Input, output, and integration a 
results are displayed in tables 
numerical integrations were done 
erenee coordinate system. 



200 

km 





400 

km 





0.0 

15 





90. 

0 rain 





30 

deg 





Noo 

n, Janua 

ry 1 , 

, 1975 



(Ful 

1 I8xl8 
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{No 

drag 
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luni-sol 
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re i 
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XIII 

(a) , 

Xlll(b) 
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XIV . 

All 

in 
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mean of 
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TABLE XIII,- TIME AND COORDINATE SYSTEM COMPARISONS 


(a) GNP versus GNP* 


Prediction interval 
davs 

, Error in position, 

km . 

GNP CPu time, 
sec 

1 .0 

negligible 

163 

2.0- 

0.00025 

328 

10,0 

0.0121 

1640 

30.0 

0.0540 

4926 

Notes : 1 . For this 

GNP and 
orbit pr 

case, the coordinate systems and time of 
GNP* coincide at the initialization of the 
ediction . 

2. In order 
the nons 

to compute the perturbing accelerations of 
pherical Earth, the satellite's postion vector 
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is rotated into an Earth-fixed coordinate system, using 
GNP or GNP». 

" i 

3. GNP* is fixed (except for the uniform rotation of Green- 
wich) whereas GNP is, changing slowly. .T.he small devi- 
ation between the two, results in the errors in .position 
shown above, 

M. The errors in position result from the effects of 
pr ece.ssion , „nuta'ti'on and time errors over 2, 

10 and 30 days. 

5. These errors remain small, growing to be only 50 meters 
at 30 days. For predictions of several months, the 
GNP* matrix should probably be updated every 3 or 
months . 


(b) GNP* versus GNP** 


Prediction interval , . Error 

days_ 

0.0623 

0.25 

0.75 

1 .0 

2.0 


in position, GNP CPu time, 

km ^ . sec 


1 .0 

■ 8' 

5.7 

30 

, 6.2 , 

88 

9.2 

117 

11 .9 

233 


Notes: 1. The GNP* matrix is computed in the same manner as 

for table XIII<a,). However, GNP** contains only 
the uniform Greenwich rotation from the mean of 1950 
coordinate system. 

2. Again, the satellite position vector is rotated from 
the inertial to an Earth-fixed coordinate system, using 
GNP* or GNP**, 

3. The position errors shown above result from the accu- 
mulated precession and nutation from 1950 to 1975 (the 
epoch of the orbit prediction,)^. For GNP**, the Earth, 
has the wrong orientation with respect to the inertial 
coordinate system. The perturbing accelerations will 
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tlierefore be 'in ‘error . The integral of these acceler- 
ation ‘ errors' results 'in the aocumulated errors sh’own 
above. 

4. The large' errors indicate that an orbit prediction pro- 
gram must account for precession-, nutation and nonuniforni 
rotation of the Earth. 

TABLE XIV.- EXECUTION TIME COMPARISONS (GNP) 


GNP GNP* GNP** 

48? ' 2 $ 0.0 


Notes: 1. Computer epu- execution times are compared here. Since 

■' GNP.** requires the least calculations, its execution 
time is taken as the reference for percentage compari-. 
sons. ' " 


2‘. The total execution times for GNP' and GNP-* are- shown 
in table X-lII{a) and Xlll(b), respectively. 

3. The GNP case requires a 4.8? increase in execution 

time whereas' GNP* requires negligible additional time. 

2,5.3 Conclusions ■- Based on these results, it Is concluded 
that the precession, nutation and nonuniform time equations can- be 
evaluated once at initialization of the orbit prediction and held 
constant thereafter’. This is the method of implementation in the 
GOPP and KSFAST programs (refs. 25 and 28), An error of about 
50 meters will accumulate (.for a nearly circular, near Earth orbit) 
after 30 days. Therefore, for- predictions over an extended time 
period, GNP* should be updated every 3 or 4 months. 


2.6 -Recommendations Based on Orbit Types 

I-t is the purpose of this section to make recommendations on 

the appropriate force models to be used for orbit prediction based 

on the type of orbit under consideration. These recommendations are 

based on the models analyzed in the previous sections. Certain 

additional force models (such as solar radiation pressure and . 

vehicle venting) may be needed f-or some satellites. 

* • < 

Ft- is difficult -t’o make -general recommendation's because of 
the large diversity of satellite 'missions' and' objectives. The 
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following recommendations are directed .primari Ijr . to the shuttle 
orbiter missions and the transfer and geosynchronous orbits of its 
payload . 

The orbit. types to be considered here are: near Earth 

orbits with small eceentrieities , elliptical transf.er orbits, and 
geosynchronous orbits. The "size and shape of each is given by 
the semimajor axis a and eccentricity e. The prediction in- 
terval to be considered is T. 


2.6.1 JJ_e.a_r-.EjLrth orbits ■- 

Orbital characteristics: 6500 km < a < 7200 km 

0 < e < 0.05 
0 < T < 2 days 

Force model recommendations for this class cf orbits: 

.a'. Geopotential - all the terms of a chosen model (such as 
ref, 42) should be- used; i.e., all the terms that were. used in the 

reduction of satellite observation data^. 

. -'b.- Atmospheric densit.y - the dynamic model AMDB* should be 
used. A small increase in accuracy may be rea-lized with the Jacohia 
model (computer cost will also be increased, however). 

c. Luni-solar gravity - these amount to about 8.0 meters in 
downrange position (after 1 day) and can usually be neglected. This 
dep.ends., .-however,, on the prediction interva.l since this error in- 
creases linearly. At- 5 days,. the error is 400 meters, 

d. Precession, ,nuta-tion, nonuniform rotation - the .GNP* 
method- should be used. 

2.6.2 . -Near-Earth orbit lifetime studies — 

Orbital characteristics: 6550 km < a < 7200 km 

0 < e < 0.05 

several- months <T<, -several y.ears 


^The order and -degree of- models can vary. It is important to 
realize, however, that each coefficient of a given model may have 
nearly equai .weight fo.r a-, near Ear th . satellite . As an- example, 
C.] 3 ^ 1 Y .may be as- physically important', as ^3 si.nce. all coeffi-, 

■ dents i = 2,3,**',N, j = 0,t,'*",i) are computed as a set. 
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Force model recommendations for this class of orbits: 

a. •Geopotential - same ^as section 2 . 6 ..1 

b. Atmospheric density - same as section 2.6.1 

c. Luni-solar gravity - Sun and Moon gravity effects should 
be included, using the analytical ephemeris model. 

d. Precession, nutation, nonuniform rotation - the GNP* 
method should be used, with update intervals of a few months. 


2.6.3 



Orbital characterlsitcs : 20 000 km < a < 25 000 km 

.65 < e < ..75 
0 < T < 5 days. 


Force mo.del recommendations for this class of orbits: 

a. Geopotential - since, the orbit .passes near the Earth, the 
complete model -should be used. Howe.ver, near apogee the higher terms 
will be insignificant and can be neglected. An automatic procedure 
for doing this has been Implemented in the program KSFAST ” (see 
ref. 28). 

b. Atmospheric density - the AMDB» or Jacchia model should 

be used, but only when. the altitude is less than about 700 km. Above 
that altitude., the drag calculation should be automatically skipped. 

c. Luni-solar gravity - the high altitude at apogee requires 
that Sun and Moon gravity be included. The analytical model of 
section 2.4 is- recommended (error is .06 km after 5 days). 

d. Precess-ion, nutation, nonunif.orm rotation - the GNP* meth- 
od should be used. 


2.6,4 


Orbital characteristics: a = 42240 km 

B =, 0-.0 

' 0 < T < 30 days 

Force model recommendations for this class of orbits: 


a. Geopotential model - The major contribution of the high 
order geopotential terms is long period resonant motion. Therefore, 
a fourth order and fourth degree model' is recommended. For pre- 
dictions of longer than 30 days,, additional tesseral terms may be_ 
needed. In any case, terms of higher order than l6-will be lost in 
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the roundoff errors of Univac 1110 double precision arithmetio, and 
should not be included. 

b. Atmospheric density - atmospheric drag can be neglected for 
this case. 

o. Luni-solar gravity - the analytical , ephemeris model may be 
used, provided that errors on the order of 0.2 km (after 30 days) 
are acceptable. ' ' 

d, 'Preoessionj nutation, nonuniform rotation - The GHP* 
model can be used-. 

2.6.5 General results .- The discussion in bee.ti,on 2.6 con- 
cerns the recommended models of atmospheric density, nonspherical 
Earth, luni-solar gravitation, coordinate systems and nonuniform 
Greenwich rotation. These force models and the three orbit types 
do not, admittedly, oover all the existing and proposed satellite 
missions-. However, the results oan be a useful help in designing 
an orbit prediction program. In the event that additional orbit 
types and force models need to be considered, a short analysis 
based on the procedures used here could produce a more accurate pro- 
gram requiring less computer runtime. Results contained in this 
report are an excellent starting point for such an analysis. 

Listed below are the important general results of these -studies 
that apply to all orbit prediction problems. 

• Periodic variations in atmospheric densi-'ty need to be in- - - 
eluded when ’air drag is important. These variations- can cause 
a 100 percent change in the effects of air drag on the orbit. 

• Numerical calculation of atmosph’erie density can be carried 
out with single precision (Univao 1110-) arithmetic. In 
addition, the solar ephemeris in a dynamic model (such as 
Jaochia and AMDB*) can be obtained from ‘a '^mean'" Sun on a 
circular orbit. The result is a 50 percent decrease in the 
computer cost of the density model evaluation, and no loss 
of accuracy. 

• The geopotential model (except for the J.g term) can also 

be evaluated in single precision without loss of accuracy. 

The geocoefficxents (table 18, ref. 42) are given to 5 
decimal digits. Therefore, intermediate calculations in 
single precision (7 to 8 decimal digits, Onivac 1110) are 

sufficiently accurate^. 

^This assumes that there are no' "programing errors" such as 
taking the sum of a large number and a small number or the dif- 
ference of two large, nearly equal numbers. ‘ 
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• The calculation of the ppe'oeasion and nutation matrices 
as well as the expression for nonuniform Greenwich rota- 
tion needs to be done only at the initialization of 

the orbit' predictions For very long predictions (several.' 
months or more) the values may need to be updated, depend- 
ing on the accuracy requirements of' the prediction* 

• For orbits passing hear the E'ar-t’h C200-400 km) it is 

necessary to' include all the known terms in the geopo- 
ten'tiai model. The reason i-s that, near the Earth, all 
the terra ‘(except hav.e about the same raagrii-tude. • 

Neglecting any term results in an error in the "mean" 
mean motion that accumulates linearly in the downrange 
direction. 

• For prediction intervals of a few days, the luni-solar 
gravitational effects ’on a near' Earth orbit (200-it,00 - km) 

•cah‘, in many cases, be' neglected . This results in typi- 
cal downrange position errors of 0,4 km after *5 days. 

With the 'shuttle orbiter, -for example, these errors are 
much less than those resulting from uncertainites in aero- 
dynamic drag. 
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